<<<<<<< HEAD ======= <<<<<<< HEAD ======= >>>>>>> f0b511cf1ec57ffdfe3f6218650b18eb2895aac2 >>>>>>> 4a484cce9de994b0c706b861618dccfbe37f8108 <<<<<<< HEAD ======= <<<<<<< HEAD >>>>>>> 4a484cce9de994b0c706b861618dccfbe37f8108 <<<<<<< HEAD ======= >>>>>>> 4a484cce9de994b0c706b861618dccfbe37f8108 <<<<<<< HEAD ======= ======= >>>>>>> f0b511cf1ec57ffdfe3f6218650b18eb2895aac2 >>>>>>> 4a484cce9de994b0c706b861618dccfbe37f8108
library(tidyverse)
library(lubridate)
library(dplyr)
library(leaflet)
library(plotly)
library(viridis)

theme_set(theme_minimal() + theme(legend.position = "bottom"))

options(
  ggplot2.continuous.colour = "viridis",
  ggplot2.continuous.fill = "viridis"
)

scale_colour_discrete = scale_color_viridis_d
scale_fill_discrete = scale_fill_viridis_d
city_files = list.files("30_cities_data")
onehundrd_city_files = list.files("100_cities_data")

Time period we are interested in

period_18 = interval(ymd("2018-02-01"), ymd("2018-04-30"))
period_19 = interval(ymd("2019-02-01"), ymd("2019-04-30"))
period_20 = interval(ymd("2020-02-01"), ymd("2020-04-30"))
period_21 = interval(ymd("2021-02-01"), ymd("2021-04-30"))

The daily mean PM2.5 AQI from Feb to Aprl of year 2019 and year 2020 in each city.

city_period_diff = 
  tibble(city = character(),
         pm25 = numeric(),
         pm10 = numeric(),
         o3 = numeric(),
         no2 = numeric(),
         so2 = numeric(),
         co = numeric(),
         )
for (city_file in city_files) {
  #print(city_file)
  
  path = str_c("30_cities_data/", city_file)
  city = strsplit(city_file, split = '-')[[1]][1]
  
  cityAir = read_csv(path) %>% 
    mutate(date = as.Date(date, "%Y/%m/%d")) %>%
    arrange(date)
  
  cityAir_19 = cityAir %>% 
    filter(date %within% period_19) 
  
  cityAir_20 = cityAir %>% 
    filter(date %within% period_20)
  
  
  pm25_19 = mean(cityAir_19$pm25, na.rm = T)
  pm25_20 = mean(cityAir_20$pm25, na.rm = T)
  pm25d = pm25_20 - pm25_19
  
  pm10_19 = mean(cityAir_19$pm10, na.rm = T)
  pm10_20 = mean(cityAir_20$pm10, na.rm = T)
  pm10d = pm10_20 - pm10_19
  
  o3_19 = mean(cityAir_19$o3, na.rm = T)
  o3_20 = mean(cityAir_20$o3, na.rm = T)
  o3d = o3_20 - o3_19
  
  no2_19 = mean(cityAir_19$no2, na.rm = T)
  no2_20 = mean(cityAir_20$no2, na.rm = T)
  no2d = no2_20 - no2_19
  
  so2_19 = mean(cityAir_19$so2, na.rm = T)
  so2_20 = mean(cityAir_20$so2, na.rm = T)
  so2d = so2_20 - so2_19
  
  co_19 = mean(cityAir_19$co, na.rm = T)
  co_20 = mean(cityAir_20$co, na.rm = T)
  cod = co_20 - co_19
  
  city_period_diff = 
    city_period_diff %>% 
    add_row(city = city, 
            pm25 = pm25d,
            pm10 = pm10d,
            o3 = o3d,
            no2 = no2d,
            so2 = so2d,
            co = cod)
}
city_period_diff = 
  city_period_diff %>% 
  mutate(
    city = paste(
      toupper(substring(city, 1, 1)), 
      substring(city, 2), 
      sep = ""),
    city25 = fct_reorder(city, pm25, .desc = T),
    city10 = fct_reorder(city, pm10, .desc = T),
    cityo3 = fct_reorder(city, o3, .desc = T),
    cityno2 = fct_reorder(city, no2, .desc = T),
    cityso2 = fct_reorder(city, so2, .desc = T),
    cityco = fct_reorder(city, co, .desc = T)
    )
city_period_diff %>% 
  plot_ly(x = ~pm25, y = ~city25, type = "bar", color = ~city25,
          colors = viridis_pal(option = "D")(3), visible = T) %>% 
  layout(title = "Feb-Aprl Daily mean PM2.5 AQI Difference, 2020 minus 2019",
<<<<<<< HEAD
=======
<<<<<<< HEAD
>>>>>>> 4a484cce9de994b0c706b861618dccfbe37f8108
         xaxis = list(title = "y"),
         yaxis = list(title = "x"),
         barmode = 'overlay',

          yaxis = list(title = "y"),
          xaxis = list(title = "x"),
          updatemenus = list(
            list(
              y = 0.8,
              buttons = list(
                list(method = "restyle",
                     args = list("x", list(~pm25)),
                     label = "pm25"),
                list(method = "restyle",
                     args = list("x", list(~no2)),
                     label = "no2")))
))
<<<<<<< HEAD
=======
======= xaxis = list(title = "PM25 AQI Difference"), yaxis = list(title = "City"))
>>>>>>> f0b511cf1ec57ffdfe3f6218650b18eb2895aac2 >>>>>>> 4a484cce9de994b0c706b861618dccfbe37f8108

Now we will see how the distribution of daily AQI differ between time period 2019 Feb-Aprl and 2020 Feb-Aprl.

city_PM25_Distribution = tibble()

for (city_file in city_files) {

  path = str_c("30_cities_data/", city_file)
  city = strsplit(city_file, split = '-')[[1]][1]
  
  cityAir = read_csv(path) %>% 
    mutate(date = as.Date(date, "%Y/%m/%d")) %>%
    arrange(date)
  
  city_19 = cityAir %>% 
    filter(date %within% period_19) %>% 
    mutate(period = "2019Feb-Aprl",
           day = format(date,"%m-%d"),
           city = city) %>% 
    relocate(city, period, day)
  
  #add a fake date "2019-02-29" with all AQI values as NA
  city_19 = 
    city_19 %>% 
    add_row(city = city, 
            period = "2019Feb-Aprl", 
            day = "02-29") %>% 
    mutate(day = as.factor(day))
  
  
  city_20 = cityAir %>% 
    filter(date %within% period_20) %>% 
    mutate(period = "2020Feb-Aprl",
           day = format(date,"%m-%d"),
           day = as.factor(day),
           city = city) %>% 
    relocate(city, period, day)
  
  city_PM25_Distribution = rbind(city_PM25_Distribution, city_19)
  city_PM25_Distribution = rbind(city_PM25_Distribution, city_20)
}
city_PM25_Distribution = 
  city_PM25_Distribution %>% 
  mutate(period = factor(period, levels = c("2020Feb-Aprl", "2019Feb-Aprl")),
         city = paste(
           toupper(substring(city, 1, 1)), 
           substring(city, 2), 
           sep = ""))
#This way is stupid but it works! I tried other methods but they just don't work as expected!!!
city_PM25_Distribution %>% 
  plot_ly(
    y = ~city, x = ~pm25, color = ~period, type = "box", 
    colors = c(rgb(0.2, 0.6, 0.8, 0.6), rgb(0.8, 0.2, 0.2, 0.6))) %>% 
  add_trace(
    x = ~pm10, visible = F) %>% 
  add_trace(
    x = ~o3, visible = F) %>% 
  add_trace(
    x = ~no2, visible = F) %>% 
  add_trace(
    x = ~so2, visible = F) %>% 
  add_trace(
    x = ~co, visible = F) %>% 
  
  layout(title = "Daily PM25 AQI Distribution, 2019 and 2020 Feb-Aprl",
         xaxis = list(title = "Daily AQI"),
         boxmode = "group",
         updatemenus = list(
            list(
              y = 1.1,
              buttons = list(
                list(label = "PM25",
                     method = "update",
                     args = list(list(visible = c(T,T, F,F, F,F, F,F, F,F, F,F)),
                                 list(title = "Daily PM25 AQI Distribution, 2019 and 2020 Feb-Aprl"))),
                list(label = "PM10",
                     method = "update",
                     args = list(list(visible = c(F,F, T,T, F,F, F,F, F,F, F,F)),
                                 list(title = "Daily PM10 AQI Distribution, 2019 and 2020 Feb-Aprl"))),
                list(label = "O3",
                     method = "update",
                     args = list(list(visible = c(F,F, F,F, T,T, F,F, F,F, F,F)),
                                 list(title = "Daily O3 AQI Distribution, 2019 and 2020 Feb-Aprl"))),
                list(label = "no2",
                     method = "update",
                     args = list(list(visible = c(F,F, F,F, F,F, T,T, F,F, F,F)),
                                 list(title = "Daily NO2 AQI Distribution, 2019 and 2020 Feb-Aprl"))),
                list(label = "so2",
                     method = "update",
                     args = list(list(visible = c(F,F, F,F, F,F, F,F, T,T, F,F)),
                                 list(title = "Daily SO2 AQI Distribution, 2019 and 2020 Feb-Aprl"))),
                list(label = "co",
                     method = "update",
                     args = list(list(visible = c(F,F, F,F, F,F, F,F, F,F, T,T)),
                                 list(title = "Daily CO AQI Distribution, 2019 and 2020 Feb-Aprl")))
                ))
            ))
## Warning: Ignoring 60 observations
## Warning: Ignoring 88 observations
## Warning: Ignoring 433 observations
## Warning: Ignoring 60 observations
## Warning: Ignoring 120 observations
## Warning: Ignoring 330 observations
## Warning: 'layout' objects don't have these attributes: 'boxmode'
## Valid attributes include:
## '_deprecated', 'activeshape', 'annotations', 'autosize', 'autotypenumbers', 'calendar', 'clickmode', 'coloraxis', 'colorscale', 'colorway', 'computed', 'datarevision', 'dragmode', 'editrevision', 'editType', 'font', 'geo', 'grid', 'height', 'hidesources', 'hoverdistance', 'hoverlabel', 'hovermode', 'images', 'legend', 'mapbox', 'margin', 'meta', 'metasrc', 'modebar', 'newshape', 'paper_bgcolor', 'plot_bgcolor', 'polar', 'scene', 'selectdirection', 'selectionrevision', 'separators', 'shapes', 'showlegend', 'sliders', 'spikedistance', 'template', 'ternary', 'title', 'transition', 'uirevision', 'uniformtext', 'updatemenus', 'width', 'xaxis', 'yaxis', 'barmode', 'bargap', 'mapType'
<<<<<<< HEAD
======= <<<<<<< HEAD
=======
>>>>>>> f0b511cf1ec57ffdfe3f6218650b18eb2895aac2 >>>>>>> 4a484cce9de994b0c706b861618dccfbe37f8108